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Abstract. In order to understand the time scales for complete mixing in a 
supernova-driven ISM we used the model of Avillez (2000) to follow the mixing 
and dispersal of inhomogeneities with different length scales. We included tracer 
fields to follow the inhomogeneities in kpc-scale simulations with different supernova 
rates. We find that an initially inhomogeneous ISM still shows local variations and 
a clumpy distribution on time scales of several hundred Myr. This clumpy structure 
does not disapear over a time of a few tens of Myr, even if the rate of SNe is 
increased to ten times the Galactic rate. These local variations occur because, even 
in the presence of numerical diffusion, gas does not mix quickly between hot and 
cold regions. The lower limit that we are able to place on the mixing time at the kpc 
scale is already longer than the time for chemistry to reach local equilibrium in the 
gas. The simulations also show that the time scale to erase inhomogeneities in the 
ISM is nearly independent of their length scale over the range 25-500 pc, contrary 
to what is expected from classical mixing length theory. 

1 Introduction 

Observations show significant abundance fluctuations along and between lines 
of sight in the interstellar medium (ISM), suggesting the ISM is far from homo- 
geneous. For example, observations using the Interstellar Medium Absorption 
Profile Spectrograph (IMAPS) show local variations in the D/H ratio of about 
a factor of three (Jenkins et al. 1999; Sonneborn et al. 2000). Deuterium 
is mainly produced during Big Bang nucleosynthesi^] and only destroyed in 
pre-main sequence stars. Therefore, the D in the ISM traces the fraction of 
gas not processed through stars, and the detection of local variations in the 
D/H ratio indicates incomplete mixing of processed material into the ISM. 
However, such an inhomogeneity picture contradicts the almost homogeneous 
distribution seen in O/H, which suggest that the ISM mixes completely on 
timescales no longer than that for chemistry reach an equilibrium. 

These observations pose the question of how chemical elements such as D 
can avoid being mixed well enough to suppress local variations in the D/H 
ratio. Furthermore, the ISM is regulated by supernova explosions that form 
well structured and explosive flows in addition to diffuse turbulence that acts 

1 However, Mullan &; Linsky (1999) have proposed that D can be formed in stellar flares 
from M dwarf stars. 



Tabic 1: Simulation parameters for the runs described in this paper. 
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a Run number and finest grid level number. 
b Finest grid resolution. 

c Time between SNe in terms of the Galactic value. 
d Time for complete mixing. 



on the smaller scales. It is therefore, necessary to understand the mixing 
process in such a medium by carrying out direct simulations of the evolution 
of the ISM polluted with inhomogencitics of different length scales. 

Section 2 deals with the description of the three-dimensional model used 
in this study, as well as, with the setup and evolution of the simulations. 
In section 3 we present the results of these simulations and discuss them in 
the light of the classical diffusion theory. Finally, in section 4 we present a 
summary of the main results. 

2 Supernova-Driven ISM Model 

2.1 Mixing Simulations 

The model includes a fixed gravitational field provided by the stars in the 
disk, and an ideal-gas equation of state. Radiative cooling is treated assuming 
collisional ionization equilibrium, using a piecewise cooling function following 
that presented in Figure 2 of Dalgarno & McCray (1972) with an ionization 
fraction of 0.1 at temperatures below 10 4 K and a temperature cutoff at 10 K. 
Background heating due to starlight varies with z as described in Wolfire et 
al. (1995), and at z = it is chosen to balance radiative cooling at 8000 K. 

Type lb, Ic, and II SNe are set up at the beginning of their Sedov phases, 
with radii determined by their progenitor masses, as described in detail by 
Avillez (2000). Sixty percent of the SNe are set up within associations, and 
the rest are set up at random sites. OB associations are set up in a layer with 



a scale height of 46 pc (from the midplanc) following the distribution of the 
molecular gas in the Galaxy, while the isolated SNe are set with an exponential 
distribution with a scale height of 90 pc. The rates of occurrence of SNe in 
the Galaxy are normalized to the volume under study, with the Galactic rate 
for types lb and Ic in the Galaxy taken to be 2 x 10~ 3 yr^ 1 , while for type II 
it is 1.2 x 1(T 2 yr" 1 (Cappellaro et al. 1997). The total rate of these SNe in 
the Galaxy is r ga i = 1.4 x 1CP 2 yr -1 , or one SN every 71 yr. 

The interstellar gas is initially distributed in a smooth disk with the vertical 
distribution of Dickey & Lockman (1990). In addition, an exponential profile 
representing the z— distribution of the warm ionized gas with a scale-height of 
1 kpc in the Galaxy as described in Reynolds (1987) is used. 

We model inhomogencitics in the ISM with a tracer field distributed in a 
checkerboard pattern of empty and full squares. We followed the evolution 
of models with squares having length L of 25, 50 and 500 pc in the x and 
y— directions (parallel to the Galactic plane). Slab symmetry of the tracer 
field is assumed along the z-direction (perpendicular to the Galactic plane). 
In the current work we report on simulations using five SN rates: a/a ga \ = 1, 
10, 20, 30 and 50. The simulations are run until 400 Myr for a/a ga \ = 1 and 
200 Myr for the remaining supernova rates. This evolution time is enough for 
the system to reach steady state. A summary of these runs is presented in 
Table 1. 

The computational domain contains a section of the Galaxy with an area of 
1 kpc 2 and vertical extension from -10 to 10 kpc. The innermost edge lies 8.5 
kpc from the Galactic centre. The computational grid has a resolution of 10 
pc, except in the layer between -250 and 250 pc, where three levels of adaptive 
mesh refinement are used, yielding a finest resolution of 1.25 pc. These simu- 
lations use the piecewise-parabolic method of Colella & Woodward (1984), a 
third-order scheme implemented in a dimensionally-split (Strange 1968) man- 
ner in combination with the AMR scheme of Berger & Colella (1989), with 
the subgridding scheme of Bell et al. (1994). 

2.2 Global Evolution of the ISM 

Once disrupted by the SN explosions, which start occurring at time t = 0, the 
disk never returns to its initial state. Instead, regardless of the initial vertical 
distribution of the disk gas, a thin disk of cold gas forms in the Galactic 
plane, and, above and below, a thick, inhomogeneous gas disk forms. The 
code does not explicitly follow ionization states, but we can trace gas with 
temperature T < 10 4 K and scale height of 180 pc, which we designate Hi, and 
gas with 10 4 < T < 10 5 K and scale height of 1 kpc, which we designate Hn. 
These distributions reproduce those described in Dickey & Lockman (1990) 
and Reynolds (1987), respectively. The upper parts of the thick Hn disk form 
the disk-halo interface, where a large scale fountain is set up by hot ionized 
gas escaping in a turbulent convective flow (see Avilez 2000 and Avillez & Mac 
Low 2001). 




Figure 1 : Time evolution of the maximum C max of the tracer field for inho- 
mogeneities with length scales of 50 pc, for a/a ga i = 1, 10, 20, 30 and 50 and 
a fine grid resolution of 5 pc. 

3 Evolution of ISM Inhomogeneities 

The tracer field is a scalar whose mass flux follows that of the density and whose 
value is unity in the regions of the flow where it is present and zero elsewhere in 
the flow. By definition the tracer field should keep its value of 1, and therefore 
its average value < C > on a checkerboard should always be 0.5 regardless of 
how it is advected. Numerical schemes introduce numerical diffusion leading 
sooner or later to a change of the value of the tracer field. However, its average 
value over the grid should remain constant. The simulations show that as the 
disk gas evolves to steady state < C > approaches 0.5, with variations of no 
more than 2.4% occurring for a supernova rate of 30 times that in the Galaxy. 
We define complete mixing to be when the maximum and minimum of the 
tracer field become numerically indistinguishable from each other and < C >, 
i.e., C max = C min =< C >~ 0.5 (AC = C max -C min < 0.01). Table 1 presents 
the time scales for complete mixing in checkerboards with 50 and 500 pc scale- 
lengths for resolutions of 5, 2.5 and 1.25 pc, while Figures 1 and 2 show the 
time variation of C max for a resolution of 5 pc and the variation of the mixing 
time with supernova rate for different finer resolutions, respectively. 

Figure 1 shows that there is a steady decrease of C max until complete 
mixing occurs. Futhermore, C max and C m i n take almost the same amount of 
time to vary from a difference AC = 1 to 0.05 as the time from AC = 0.05 
to 0.01, when we define complete mixing to occur. This indicates that smaller 
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Figure 2: Variation of the mixing time with supernova rate for finer grid 
resolutions Ax — 1.25, 2.5, 5 pc. Mixing time increases by factor between 
1.05 and 1.2 when resolution is incresed by a fetor of 2. 



scale structures dominate the mixing process after the gas has been completely 
mixed on larger scales. 

Figures 1 and 2 show that the time taken for complete mixing to occur 
decreases with increasing supernova rate, that is, increased energy input into 
the ISM. However, the rate of increase of mixing reaches a saturation level 
for oloG a i axy > 20, beyond which the increase in the supernova rate does not 
accelerate mixing in the ISM. For SN rates greater than 20 times the Galactic 
rate we see that the mixing time scales differ by only a few Myr (see Table 
1), with complete mixing occurring at ~ 107 Myr for Ax = 5 pc. This results 
from the fact that for a/acaiaxy > 20, after the first 50 Myr of evolution a 
large fraction fo the disk gas has a temperature greater than 10 5 K and at 100 
Myr of evolution all the gas is hotter than 10 5 K. As a consequence of the high 
SN rate the disk becomes hotter and cooling becomes inefficient. At this high 
temperature the sound speed is correspondingly high, so the gas mixes within 
a 40-50 Myr period (Roy & Kunth 1995) leading to complete mixing at some 
110 Myr. 

The time taken for complete mixing to occur in a checkerboard with 50 
pc squares varies from ~ 300 Myr for u/u Galaxy = 1 to some 120 Myr 
for o I ctGaiaxy = 20 and 107 for u/ucaiaxy > 30 for Ax = 5 pc, while for 
Ax = 1.25 pc mixing varies from 347 Myr for cr/acalaxy = 1 to 150 Myr for 
&l&Gaiaxy = 20. On 500 pc squares complete mixing is surprisingly only de- 
layed by around 20 Myr (Table 1). The similarity between these time scales 



indicates that the mixing process is almost independent of the length-scale of 
the inhomogeneities, and therefore, of the initial setup adopted for the tracer 
field. Rather, it is regulated by supernovae and their energy input into the 
ISM. 

According to classical mixing length theory the mixing time scales of these 
boards should follow the relation r cx £ 2 , where I is the length of the squares. 
Therefore, the mixing time of inhomogeneities with a scale length of 500 pc 
should be 25 times longer than the mixing time for 50 pc squares. The com- 
puted similarity in the mixing time scales (Table 1) for the different inho- 
mogeneity sizes shows that classical mixing length theory does not describe 
mixing in the ISM well. This seems to be because much of the ISM is dom- 
inated by nearly laminar flows inside of SNRs and superbubbles rather than 
by diffuse turbulence. The latter does act on smaller scales and in cooler gas, 
in which classical mixing length theory may apply better. 

4 Summary 

In this paper we try to shed light on the time scales that it takes to erase 
inhomogeneities in the ISM and how this process depends on the inhomogene- 
ity length scale. We find: (1) The time scale to erase inhomogeneities in the 
ISM is nearly independent of their length scale on scales above 10 pc or so; 
(2) Mixing time scales decrease with increasing supernova rate until all the 
gas is hot and the mixing rate saturates. This occurs when c/cQalaxy ^ 20, 
when the mixing time scale is ~ 107 — 120 Myr for Ax = 5 pc and 150 Myr 
for Ax = 1.25 pc; (3) Even after some 80% of the gas has been mixed, the 
simulations show that the ISM remains clumpy for at least twice as long if we 
examine scales smaller than a kpc. Even if the rate of SNe is increased to ten 
times the Galactic rate, this clumpy structure does not disappear over a time 
of a few tens of Myr. 

This work shows that inhomogeneities that are present in the ISM, assum- 
ing that no further inhomogeneities are introduced into the system, take a 
very long time to be erased. The time for complete mixing (greater than 107 
Myr, which occurs for cr/o'Q a \ aiX y = 50 and Ax = 5 pc) is already longer that 
any time scale for chemical equilibrium to occur, and therefore, the ISM will 
show a clumpy distribution. The local variations occur because, even in the 
presence of numerical diffusion, gas does not mix quickly between hot and cold 
regions. 
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